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Abstract 


We  study  a  two-dimensional  reconstruction  algorithm  due  to  D.C.  Barber  and 
B.H.  Brown,  applied  to  a  linearized  electrostatic  inverse  problem.  Firstly  we  demon¬ 
strate  how  this  algorithm  fits  within  the  framework  of  inverses  of  generalized  Radon 
transforms  studied  by  G.  Beylkin.  Secondly  we  construct  an  iterative  improvement 
of  the  Barber-Brown  algorithm  based  on  the  conjugate  residual  method.  We  present 
several  numerical  results  obtained  with  this  iterative  algorithm _ _  r'' 


J 


1  Introduction 


In  electrical  impedance  imaging  one  seeks  to  reconstruct  the  internal  conductivity  (or 
impedance)  profile  of  an  object  from  boundary  measurements  of  voltages  and  corre¬ 
sponding  current  fluxes.  There  has  been  significant  advances  made  in  recent  years  on 
both  practical  aspects  of  the  reconstruction  problem  [2,  3,  11,  16,  17]  as  well  as  on  the¬ 
oretical  aspects  of  the  uniqueness  and  continuous  dependence  question  [1,  10,  14,  15, 

18].  We  shall  not  attempt  a  review  of  the  literature;  the  reader  may  consult  [6]  for  an 
extensive  list  of  references. 

One  particular  algorithm  which  haw  shown  itself  to  be  surprisingly  effective  given  its  i 

low  cost  was  developed  by  D.C.  Barber  and  B.H.  Brown  [2,  3]  for  use  in  the  context  of 
medical  tomography.  In  this  paper  we  analyze  the  Barber  and  Brown  method  in  some  | 

detail.  In  particular  we  show  that  the  crucial  backprojection  component  of  the  algorithm 
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may  simply  be  viewed  as  part  of  an  approximate  inverse  of  a  generalized  Radon  transform, 
as  studied  extensively  by  Beylkin  [4,  5], 

The  Barber  and  Brown  algorithm  has  already  been  compared  numerically  to  several 
other  algorithms  [19]  and  appears  quite  exceptional  in  achieving  a  moderate  accuracy  at 
an  extremely  low  cost.  In  the  latter  part  of  this  paper  we  construct  an  iterative  extension 
of  Barber  and  Brown’s  algorithm  based  on  a  conjugate  residual  method.  This  new  algo¬ 
rithm  can  achieve  a  significantly  higher  accuracy  at  a  modest  increase  in  computational 
cost. 

In  the  analysis  and  computations  of  this  paper,  we  restrict  ourselves  to  a  circular 
domain  and  consider  only  the  linearized  identification  problem.  We  are  currently  in  the 
process  of  implementing  the  iterative  algorithm  for  the  full  nonlinear  problem  on  an 
arbitrary  polygonal  domain.  Results  from  this  work  will  be  reported  elsewhere. 

The  paper  is  organized  into  seven  sections.  In  Section  2,  we  give  a  brief  description 
of  the  electrical  impedance  problem,  and  review  Barber  and  Brown’s  backprojection 
scheme.  In  Section  3  we  show  how  the  Barber  and  Brown  backprojection  fits  within 
the  framework  of  approximate  inverses  of  generalized  Radon  transforms,  constructed  by 
Beylkin.  We  also  examine  the  relation  between  the  linearized  elliptic  forward  problem  and 
a  generalized  Radon  transform.  Section  4  provides  a  discussion  of  our  implementation  of 
the  Barber-Brown  algorithm.  In  Section  5  we  study  properties  of  the  ronductivity-to-data 
map  and  its  composition  with  the  backprojection;  in  particular,  we  perform  an  eigenvalue 
analysis  to  study  the  sensitivity  of  these  maps.  This  study  leads  to  the  development  of 
an  iterative  algorithm  based  on  a  conjugate  residual  method,  as  presented  in  Section  6. 
Results  from  our  numerical  experimentations  with  the  iterative  method  are  displayed  in 
Section  7. 


2  The  linearized  inverse  problem-A  backprojection 
reconstruction 

As  a  mathematical  model  for  the  direct  current  electrostatic  problem,  we  suppose 


V  •  (7VU)  =  0  in  fi  (1) 

du 

y-r-  =  tf>  on  cm, 
an 

where  u  is  the  voltage  potential,  ^  is  a  boundary  current,  and  y  denotes  the  conductivity 
profile.  The  linearized  problem  for  a  small  perturbation  S y  (of  7)  and  corresponding 
perturbation  6U  (of  U)  now  becomes 


V  •  (7V6C/)  =  -V  •  (6yVU)  in  fl, 

d(6U)  dU 

r—  =  on  da 

In  the  remainder  of  this  paper,  we  shall  for  simplicity  assume  that 


(i)  Q  is  the  unit  ball,  0  =  {a:  G  IR2  :  |x|  <  1}. 

(ii)  7=1. 


slon  For 

GRAM- 
TAB 

Minced 


□  □ 


where  dbu/dr  is  the  counter  clockwise  tangential  derivative  of  a  Dirac  delta-function, 
bu,  u  6  dQ.  We  furthermore  assume  that 


(iii)  b 7  =  0  near  the  dipole. 


Under  assumptions  (ii)  and  (iii)  the  linearized  problem  reduces  to 


A  bU 
d(6U) 
dn 


-V(by)-VU  in  fi, 
0  on  dSl.- 


(3) 


The  original  problem  in  electrical  impedance  imaging  is  to  determine  a  consistent  y, 
given  knowledge  of  u|an  (the  solution  to  (1))  for  various  choices  of  t p.  The  linearized 
inverse  problem  associated  with  (3)  therefore  becomes 

(I  P)  Given  bU\dn  for  various  choices  of  dipole  solution  U ,  determine  a  consistent  in¬ 
crement  by. 


We  use  the  notion  consistent,  since  of  course  none  of  these  inverse  problems  can 
be  expected  to  have  a  unique  solution  given  only  a  limited  set  of  choices  of  boundary 
fluxes.  The  algorithm  of  Barber  and  Brown  represents  an  ingenious  one  step  approximate 
solution  to  the  problem  (I  P).  To  explain  the  algorithm  we  briefly  review  its  derivation 
(along  the  lines  of  [3]). 

Since  f2  is  the  unit  ball,  the  solution  to  (2)  is  known  in  closed  form: 
x' 

u  =  -irrz 7S  with  *i  =  ■  *>  *2  =  1  -  w  ■  *. 

Xl  T  X2 

where  w1  =  (— W2,  uq)  is  the  rr/2  rotate  of  the  dipole  location  w  =  (wi,  The  function 

V  =  „  x'2 

is  a  harmonic  conjugate  to  U  on  f2.  Indeed,  x  — ►  (<7,  V)  (the  so-called  Poincare  map) 
conformally  maps  D  onto  the  upper  half  plane  P  =  {V  >  1/2}.  The  problem  (3)  simplifies 
in  the  (U,  V)  coordinates  to  read 


A  6U 

d(6U) 

dV 


d(by) 

dU 


in  P 


0  on  dP  =  {V  =  1/2}. 


(4) 


Note  that  the  function  by  =  by(x(U,V,u))  is  now  a  function  of  U,  V  and  w.  The  extra 
data  being  used  to  reconstruct  by  is  the  function  bU |v  =  i/2  (a  function  of  U  and  u).  For 
a  single  fixed  dipole  location  u>o,  a  consistent  conductivity  increment  by0  is  given  by 

ho(x)  =  -Qy(6V\v_i/2)(U(x,uo),uo), 


as  follows  easily  from  (4).  Barber  and  Brown  suggest  the  average 


\  f  & 

B(x)  =  ~2n  !  i(^/<5t/lv'  =  i/2)(s-w)l»=t/(r,a,)<J>(x,w)d5u 


(5) 


a  rough  approximation  to  the  conductivity  increment  by  ( D  is  oi  course  not  in  gciieial 
consistent  with  any  dipole  measurements). 


The  formula  (5)  has  a  geometrical  interpretation  illustrated  by  Figure  1.  Given  a 
point  x  to  be  imaged,  and  a  dipole  location  w,  consider  the  equipotential  circular  arc 
{z  :  U(z,ui)  =  U(x,u)  =  s}  which  originates  at  u>  and  passes  through  x.  The  point 
where  this  equipotential  arc  intersects  the  boundary  is  z(s,  1/2,  w)  (in  U,  V  coordinates, 
for  fixed  w,  it  corresponds  to  U  =  s,  V  =  1/2).  The  first  term  in  the  integrand  is  the 
known  quantity  (-§?&U/ £fU)(x(s,  l/2,w),u>)  (see  (6)  below). 


U(z,u)  =  U(x,u)  =  s} 


Figure  1  In  this  picture,  the  dipole  is  located  at  ui,  the  equipotential  arc  through  x  is  displayed. 
The  intersection  of  the  equipotential  arc  with  the  boundary  is  r(s,  1/2,  u>). 


The  weight  $  is  selected  in  a  very  particular  fashion.  Let  p  be  the  angle  between  the 
tangent  to  the  curve  U(z,w)  =  U(x,w)  at  x  and  the  vector  (0, 1).  Denote  by  9  the  angle 
between  (0,-1)  and  w  (cf.  Figure  1).  The  function  $(x,w)  equals  \8p/89(x,  w)|,  so  that 

*(x,U)dSlAl  =  \^(x,u,)\d9  =  dp, 

(for  a  fixed  imaging  point  x).  Therefore,  the  average  in  (5)  is  exactly 

r2r 


~2n  Jo  ^^7l5t/lv,=1/2>(s’u') 


which  corresponds  to  a  uniform  distribution  of  angle  p  (not  9). 

It  is  not  difficult  to  see  that 

$(s,w)  =  2V(x,  w)  -  1 

and 

(^«/k=i/2)(^)  =  (±6U/^U)(x(s,  1/2,  w),  w). 

Indeed,  in  Figure  1,  p  —  p'  —  9  and  so 

dp  dp' 

89  ~  89 

On  the  other  hand  in  the  x'  coordinates  (for  fixed  w)  the  tangent  to  the  curve 
{z  :  U(z,u )  =  U(x,u)}  is 


(ба) 

(бб) 


aud  therefore 


-'2  _  x '2 
I  *2  *  1 

coap  =  WT*F 


sin  p 


2z',zi 

zF+zF' 


(7) 


Differentiation  of  the  first  formula  in  (7)  yields 

d£__d_  +  »? 

30  30  \  x'2  +  x'2  )  dO  Vzf  +  /  2*1X2 


and  a  simple  computation,  using 


now  gives 


Consequently 


d  ,  _  ,  ,  d  ,  _  , 

^2-  *1. 


-  =  2-2V. 

60 


which  immediately  leads  to  the  formula  (6a). 

The  formula  (6b)  represents  a  simple  change  of  variable,  since 

(£j6U\v=1/2)(s,u)  =  XcT(z(*, 1/2- w)'w)  =  1/2,  W),  W), 

where  3/3r  is  the  counter  clockwise  tangential  derivative  and  x(s,  1/2,  u)  is  as  in  Figure 

1. 


The  Barber  and  Brown  average  may  now  be  written 

B(x)  =  ±.  J  {±su/£u)(x(s,  l/2,U),w)\tmUl,„)(l-2V{z,U))dSu.  (8) 

Except  for  a  filter,  which  we  shall  briefly  discuss  later,  the  calculation  of  the  integrals 
(8)  represent  exactly  the  reconstruction  method  suggested  by  Barber  and  Brown.  In  the 
following  section,  we  explain  how  (8)  may  be  seen  as  part  of  an  approximate  inverse  for 
a  generalized  Radon  transform. 


3  The  generalized  Radon  transform  and  its  inverse 

The  function 

*(*.0  =  lf|tf(*.€/KI).  *en,  ZeiR2\{0}. 

is  positive-homogeneous  of  degree  one: 

*(*,  A0  =  A*(*,0,  (9) 

it  is  also  infinitely  often  differentiable  in  D  x  (ZR2  \  {0})  and 

v,*(*,o#o  v*efU€ZR2\{ 0}.  (io) 

The  function  <j>  defines  a  family  of  arcs  (parts  of  circles)  to  be  used  for  the  generalized 
Radon  transform 


H,tUl  =  {z  €  D  :  <t>{x,u)  =  s},  s  €  ZR,  |w|  =  1. 


5 


As  a  measure  ou  each  arc  H,  w  we  take 


dfi  =  \Vt<t>\~ld<r, 

where  dtr  denotes  standard  arclength;  we  let  a  denote  the  amplitude 

a(x,u))  =  IV^x.w)!2. 

Following  the  notation  in  [4]  we  now  define  a  generalized  Radon  transform 
(f?u)(s,w)  =  /  u(x)a{x,ui)dp 


(11) 


(12) 


■  L 


u(x)\V  X<j>(x,u>)\dcr 


for  any  u  €  Co°(fl). 

Beylkin  provides  a  recipe  for  the  approximate  inversion  of  transforms  like  that  in 
(12).  His  main  assertion  is  that 

R*KR  =  id  +  T,  (13) 

where  T  is  compact:  L2(f2,  compact)  — »  L2(fl,  loc)  and  R *  is  the  so-called  backprojection 


(**«)(«)=  [  ^\v(s,u >)\t=4lx<u)dSu. 

•/  |u/|  =  1  ®(*tw) 


K  consists  in  convolution  with  the  generalized  kernel 


The  function  h  is  given  by 


=  2^ /I1'1™' 


(14) 


(15) 


The  weight  h(x,  v)/a(x,  w),  |w|  =  1,  appearing  in  R*  has  a  very  simple  relation  to 
\dp/dO\  =  2V  —  1  (cf.  Figure  1). 

Lemma  1 


h(x,u)/a(x,w)  =  2V(x,w)  —  1,  x  €  fi,  |w|  =  1. 


Proof 

Using  polar  coordinates  (r^,0^)  in  the  {-plane,  we  get 

h(x,u>)  =  [(i-V^^fiv,#,.) 

for  |w|  =  1.  Since  0  is  positive-homogeneous  of  degree  1  and  since  6 j  is  the  same  as  the 
angle  9  (corresponding  tow  =  {/|{|  in  Figure  1)  we  get 

A(x,w)  =  [(Vr^-L(Av^)](x,w).  (16) 

Note  that  the  formula  (16)  only  involves  the  function  <f>(x,w),  |w|  =  1,  due  to  the  fact 
that  differentiations  with  respect  to  6  and  x  are  performed  with  r ^  =  |{|  =  constant  =  1. 


As  before,  let  x'  denote  the  coordinates  related  to  x  by  the  orthogonal  affine  trans¬ 
formation 

-'-«-+(?)•  zi) 

Then  Vr  =  Qj^X‘,  and  from  (16)  we  therefore  get 


h(x,u>)  =  (QjVr-t)1  •  ((^g»)TVt^)  +  (V,^  •  (§gVr‘4>) 

=  |v^l2  +  (vr«0)M^v^). 


Division  by  a  =  |Vr<£|2  =  |Vr<<£|2  yields 


h(x,u)) 

a(x,u) 


=  1  + 


(Vr^)x  (&VX><*) 


I  vr^|  |Vr^| 


(2,w). 


(17) 


(18) 


The  function  <j>  satisfies 


*<*■“>  =  ■ 


and  dx[/dO  —  x'2  —  l^dx^/dO  =  — x(.  A  simple  calculation  based  on  (18)  now  leads  to 


h(x,u)/a(x,u)  =  2— 


x2 


x'l  +  *2^ 


-  1 


=  2V(*,w)-l. 


□ 


Remark  3.1 

We  note  that  in  his  paper  [4],  Beylkin  assumes  that  the  phase  function  <j>  is  odd 
with  respect  to  the  variable  £  (he  assumes  that  <t>  is  not  just  positive-homogeneous,  but 
homogeneous  of  degree  one).  The  only  place  this  is  used  in  an  essential  way  is  in  the 
calculation  of  the  splitting  of  the  Fourier  Integral  Operator  F,  on  page  583.  F  is  given 
by 


with 


and 


Fv(y)=  f  G(y,u)dSu , 

-f|w|  =  l 

G(y,u)=  — ^  J  e,*^x'y,rw^A(x,  y,  ru)u(x)dx^j  rdr, 


a(x,  rtr)h(y,£) 

^(*.y.O  =  4>{x,0  -  <Kv,t),  A(x,y,o  =  — r— n — 

a<*  ifr> 

If  4>  is  not  necessarily  odd  with  respect  to  £,  then  a  slightly  different  calculation  gives 

and  therefore 


provided  ti  is  real.  In  other  words 


MV 

V 


-Si 


‘to 

•M 


* :»**:■ 


.‘V 


•v 

*.*!• 


:;v.a 


Re(G(y,w))  =  y^KR(umy,U)), 
a(y.w) 

where  K  represents  convolution  with  the  generalized  kernel 


—  fl 
2*)2  J- oc  ' 


r|e,r’dr. 


The  formula  that  corresponds  to  (3.3)  on  page  584  of  [4]  in  this  case  becomes 

Re(Fxt(y))  =  R*I<R(u)(y),  (19) 

provided  ti  is  real.  There  are  no  changes  required  to  show  that 

F  =  id  +  T  (20) 

where  T  is  compact:  L2(fi,  comp )  — ►  L2(Q,  loc).  Based  on  (19)  and  (20)  we  immediately 
conclude  that 

R*KR  =  id+T 

where  T  is  compact:  L2(Q,  comp) —*  L2(Q,  loc).  □ 

Let  W  denote  the  function 

From  Lemma  1  we  immediately  get 
(**W)(x) 

=  /  M*>w)/a(*-w)W(s,w)|J=*(t|W)dSu, 

=  ^r/w|_1(|:w/|:y)(x(s’i/2-u,)>u')l*=^.‘-)(i  -2V(x,u,))dSu. 

We  shall  shortly  verify  that 

W  ss  KR(6j)  (  in  a  very  crude  sense).  (21) 

We  therefore  conclude  that  the  Barber  and  Brown  average 

^/«l-i(^/^)(r(S’1/2’W)’W)!'=i;(t,,,')(1"2K(l’a,))d5u'  =  R'W{X) 

ss  R*KR(Sy)(x) 

ss  6y(x) 

may  be  viewed  as  a  crude  first  approximation  to  67(x).  We  find  it  quite  remarkable  that 
Barber  and  Brown  by  purely  heuristic  arguments  were  able  to  find  the  “right”  weight 
{2V  —  1)  for  their  backprojection. 


We  complete  this  section  by  showing  (21).  The  function 


G(u°v°>iu'V)  =  -b{tu^u-! 


U-U0 


U-Uo 


)2  +  (v-  vQy  ^  (u  -  u0y  +  (v  +  v0-i y-)' 


Vo  >  1/2,  solves 


&G(Ua,va)  =  -Qfj^Uo.Vo)  *n  F=  Rx  (l/2,oo) 

dG'{dV~  =  0  on  dP={V=  1/2}. 


It  follows  from  this  and  (4)  that 


6U{U,V,u)=  [°°  [°°  G(U,tV,)(U,V)6y(U',V',cj)dV'dU'l 
J-oo  J 1/2 


and  consequently 


-(~SU\v=l/2)(s,u) 


=  -  j  J  ^jG(Wiv){s,ll2)6-i{U'  X^dV'dU' 


Fourier  transformation  of  (22)  with  respect  to  s  leads  to 


-{■£j6U\v=1/2)A(r,u>)  =  I|r|  ^“e-(v'-1/2)l’-|^A(r,K',W)dK/; 


here  we  have  used  the  fact  that 


((&*)'<■>-/. 


»  a*  -  U 2 


dU  =  Trlrle" 


q  >  0. 


V(«3  +  ^a)V  v  /  7-oc  (g2  +  t/2)2  - -  ’  " 

If  we  replace  e-^  -1/2)lrl  by  its  value  at  r  =  0  (or  V'  —  1/2)  then  (23)  reads 


or  equivalently 


-(^tf|v=i/2)A(r,w)  *  i|r|  (jH  tf7(-,  K',W)dl/')A(r), 

-  (^I7|k«i/3)(*,«)  *  2*rA:(^~ 


where  K  represents  convolution  with  the  generalized  kernel 


—  r 

2*)2  L C 


|r|e,r’dr. 


According  to  (6)  the  left  hand  side  in  (24)  is  exactly 

~(§;8U/ §^U){x(s,  l/2,w),w), 

and  by  a  simple  change  of  variables 

/  «7(i,K'lW)dF'=  /  «7(x)|V,Cf(x,w)|d<r  =  fl(«7)(«,( 

>'1/2 


Therefore 


as  stated  in  (21). 
Remark  3.2 


ss  A' A(57)(s,w) 


We  note  that  the  approximation  (21)  is  best  for  smooth  6j  or  6y  whose  “singularities” 
lie  not  too  far  from  the  boundary.  For  6-y  with  “singularities”  near  the  center  (21)  repre¬ 
sents  only  a  very  crude  approximation.  In  contrast  the  approximation  R" K R(6y)  ss  <*> 7 
seeks  to  fit  “singularities”  and  represents  only  a  very  crude  approximation  for  smooth 
6y.  This  difference  in  the  nature  of  the  approximations,  we  believe,  is  the  main  reason 
that  the  Barber-Brown  backprojection,  R*W ,  in  itself  only  gives  a  crude  approximation 
to  fi7.  □ 


4  Implementation  of  the  backprojection  and 
the  Barber-Brown  filter 

This  section  contains  a  description  of  the  numerical  implementation  of  the  backprojection 
discussed  in  the  previous  section,  including  a  brief  description  of  the  spatially  dependent 
filter,  which  is  used  to  improve  the  reconstruction. 

Let  xb  denote  the  point  xh  —  z(s,  1/2,  w)  on  the  boundary  of  the  unit  circle  obtained 
by  solving 

U(x,ui)  =  s,  with  |x|  =  1,  (25) 

for  a  specified  s  and  ui.  The  data  based  on  which  we  seek  to  reconstruct  67  is 

W(s,u)  :=  -lsu/lu(xitu). 

The  backprojection  formula  (8)  amounts  to 

B(x)  =1  [  W(s,c  )\,^xm)(2V{xm)  -  1  )dSu. 

The  discrete  approximation  to  the  back  projection  is  represented  by  a  matrix,  for 
simplicity  also  denoted  B.  For  the  iterative  scheme  we  will  need  the  entire  matrix  B  not 
just  its  action  on  individual  vectors. 

We  assume  that  the  experimental  setup  contains  m  electrodes.  The  midpoints  be¬ 
tween  electrodes  are  numbered  1  through  m  (see  Figure  2).  A  pair  of  adjacent  electrodes 
through  which  input  current  flows  (in  and  out  respectively)  is  called  the  driver  pair.  As 
a  realization  of  a  dipole  at  location  1,  we  select  the  driver  pair  to  be  the  two  electrodes 
adjacent  to  location  1.  We  measure  voltage  differences  on  all  electrode  adjacent  pairs. 
The  rescaled  values  of  these  differences  represent  our  discrete  data  corresponding  to  a 
dipole  at  location  1.  As  pointed  out  by  Barber  and  Brown  it  is  often  not  possible  to 
obtain  reliable  values  of  the  voltage  tangential  derivatives  at  the  dipole  location,  and  the 
locations  adjacent  to  it  (in  this  illustration,  locations  1,  2  and  m).  This  is  because  mea¬ 
surements  of  voltages  at  the  driver  pair  tend  to  be  inaccurate.  They  suggest  filling  this 
data  gap  by  extrapolation  from  neighboring  measurements.  In  our  numerical  simulation, 
we  assume  that  the  tangential  derivatives  of  the  voltage  have  been  obtained  at  locations 
adjacent  to  the  dipole  location,  and  use  the  fact  that  the  quantity  \V  in  the  limit  of  a 
perfect  dipole  is  zero  at  the  dipole  location  (see  (22)).  To  obtain  the  complete  data  set, 
we  cycle  through  all  possible  driver  pairs. 

The  integral  above  is  replaced  by  a  sum  over  all  driver  pairs.  Thus  the  discrete  version 
of  the  backprojection  is 

1  W^s^i)\>=V{x^,){^V{x,u)j)  -  1).  (26) 

j  —  1 ,  m 

Notice  that  the  point  xb  on  the  boundary,  which  satisfy  (25)  for  s  =  U(x,ujj)  may  not 
lie  at  a  measurement  point.  Hence  interpolation  of  the  data  may  be  necessary. 


Figure  2  The  pixel  being  illuminated  has  index  t  with  corresponding  location  x.  The  driver 
electrodes  are  the  pair  adjacent  to  1.  This  experiment  generates  a  dipole  at  1.  Collected 
measurements  at  the  numbered  locations  produce  the  data  { w\ ,  w2,  W3,  ■  ■  ■ }  dipoic  at  1  •  The 
equipotential  through  x  intersects  the  boundary  at  xb  at  angle  8b.  The  indices  of  the 
electrodes  nearest  xb  are  l  and  1  +  1. 

To  represent  the  data  discretely,  we  use  an  m2  vector  w.  The  index  on  w  locates  a 
dipole  and  a  point  on  the  circle  at  which  the  tangential  derivative  of  the  potential  is 
measured.  The  value  at  that  index  represents  the  quantity  W .  The  vector  w  consists  of 
m  m-vectors.  Each  vector  corresponds  to  one  experiment,  that  is,  it  contains  the  ratios 
between  the  measured  tangential  derivatives  and  reference  tangential  derivatives  at  all 
measurement  locations  for  a  fixed  dipole.  Thus 

W  xz  [(wi ,  W21  U^3,  .  .  ■ )dipole  at  1 ,  (u>l ,  W2^  U>3,  •  ■  -)dtpole  at  2 1  ( •  •  •) dipole  at  3i  '  '  ']  • 

The  original  circular  domain,  |z|  <  1,  is  embedded  into  [— 1, 1]  x  [—1, 1]  for  convenience 
of  graphics.  The  extended  domain,  which  we  call  the  image  domain,  is  discretized  so  that 
there  are  n2  pixels  of  size  2/n  by  2/n.  The  conductivity  perturbation  6y(x)  is  replaced 
by  an  n2  vector,  for  simplicity  also  denoted  by  6y.  The  index  on  the  vector  6y  locates 
a  pixel  in  the  image  domain,  and  the  value  at  that  index  represents  the  conductivity 
increment  at  the  center  of  the  pixel. 

-  2  2 

The  backprojection  matrix  B  takes  a  vector  in  IRm  and  maps  it  to  JR"  .  The  row 

vectors  of  B  which  correspond  to  pixels  outside  of  |z|  <  1  are  set  to  zero.  We  compute 
the  matrix  B  by  rows.  To  describe  the  computation  of  B,  let  us  take  a  pixel  indexed  by 
»,  centered  at  x.  The  procedure  outlined  below  is  based  on  formula  (26). 

For  illustrative  purposes,  let  us  take  the  case  when  the  dipole  is  located  at  index  j  —  1. 
With  this  dipole  location,  we  can  calculate  the  first  ro-entries  of  the  ith  row  of  B.  Our 
computational  experience,  as  well  as  that  of  Barber  and  Brown,  suggests  that  a  certain 
smoothing  is  required.  This  smoothing  is  accomplished  by  imaging  several  points  xk , 
neighboring  z,  and  then  replacing  the  terms  W{s,Wj)\)=U(z,w,)  in  (26)  by  the  respective 
averages  over  the  points  xk .  As  a  common  weight  we  use  the  value  of  ( 2V  -  1  )/m  at  the 
point  z. 
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The  weight  at  x  is 


wt  = 


1,  2(*2  +  l) 

2iryx\  +  (z2  +  l)2 


A  quick  calculation  shows  that  the  equipotential  arc  connecting  the  dipole  location  1  and 
xk  intersects  the  circular  boundary  at 


From  this,  we  calculate  the 


x\  =  4Uk  /(4(Uk)2  +  1), 

=  2/(4(Uk)2  +  1)  —  1. 
angle  6b  to  be 

8b  —  tan~l{—x\/ x  2). 


It  is  now  simple  to  identify  the  indices  /  and  /  +  1  corresponding  to  measurement 
locations  between  which  the  point  xb  lies.  Linear  interpolation  gives  factors 

fi  =  (Oi-0,)/ A0, 

/1+1  =  (®i+i  —  0i)/A0, 
corresponding  to  /  and  /  +  1  respectively. 

In  summary  our  algorithm  for  constructing  the  first  m  entries  of  the  »th  row  of  the 
matrix  B  is: 

1)  for  each  point  xk 

a.  initialize  the  m-vector  workk  to  zero 

b.  compute  6b  and  identify  indices  /  and  /  +  1 

c.  in  the  vector  workk  set  the  1th  and  the  /+lth  entries  to  //  and  /(+1  respectively 

2)  form  the  average  of  all  the  vectors  workk  and  store  result  in  the  vector  work 

3)  calculate  the  weight  wt  and  multiply  work  by  wt 

The  resulting  vector  comprises  the  first  m  entries  of  the  *th  row  of  B.  The  scalar  multipli¬ 
cation  of  this  row  vector  with  the  first  m  entries  of  a  data  vector  w  is  the  contribution  of 
the  (averaged)  operation  (26)  for  j  =  1  at  the  image  point  i.  To  find  the  next  m  entries, 
move  the  dipole  to  location  2  and  carry  out  the  same  calculation. 

In  our  implementation,  we  take  m  =  n  =  16.  However,  in  order  to  obtain  an  accurate 
B,  we  initially  calculate  it  for  m  =  6A  The  resulting  matrix  is  postmultiplied  with  an 
interpolation  matrix  such  as  to  make  B  a  256  x  256  matrix.  The  interpolation  is  linear. 
To  achieve  a  desirable  level  of  smoothing,  we  take  25  points  xk ;  these  points  make  up  a 
regular  stencil  covering  an  area  the  size  of  the  pixel  centered  at  x. 

The  matrix  B  was  used  in  Barber  and  Brown  [3]  to  reconstruct  point  images  at  various 
positions  in  the  circle.  It  was  noted  that  the  reconstruction  had  limited  resolution,  and  it 
was  also  noted  that  the  resolution  depended  on  the  position.  In  order  to  focus  the  recov¬ 
ered  image,  they  designed  a  position  dependent  filter.  The  construction  of  this  filter  is 
purely  heuristic,  and  some  of  the  parameters  are  arrived  at  by  experimentation.  However, 
we  cannot  overlook  its  effectiveness  and  so  we  have  included  it  in  our  implementation. 


Briefly  described  the  filter  works  as  follows:  for  a  fixed  pixel  at  location  x  one  con¬ 
structs  a  combination  of  a  rotation  and  a  conformal  map  in  order  to  map  the  unit  ball 
to  itself  and  map  x  to  the  origin  (the  rotation  takes  x  to  |x|  on  the  positive  real  axis, 
the  conformal  map  is  now  z  — <•  (2  -  |x|)/(l  -  z  |x|),  using  the  natural  identification  of 
Bi2  and  (T  ).  In  the  new  copy  of  the  unit  ball  one  computes  the  effect  of  convolution 
of  the  transformed  Btu-values  with  a  Gaussian  distribution  centered  at  the  origin  (one 
actually  only  computes  the  convolution  at  the  origin,  the  point  which  corresponds  to 
location  x).  Finally  the  filtered  Bta-value  at  the  location  x  is  taken  to  be  c\  (original 
Bt/>-value)— C2  convolution.  The  constants  c\  and  c 2  are  empirically  chosen  to  be  16  and 
15  respectively!  For  more  details  we  refer  the  reader  to  [3], 

We  used  Barber  and  Brown’s  code  without  changes  to  calculate  this  filter  matrix, 
which  we  denote  by  F.  The  premultiplication  of  B  by  F  gives  what  we  call  the  filtered 
Barber-Brown  backprojection,  denoted  here  by 

B  =  FB. 

Thus  given  a  data  vector  w,  corresponding  to  an  experiment,  we  find  a  rough  reconstruc¬ 
tion  through  67  ss  Bw.  Numerous  simulated  reconstructions  using  the  matrix  B  can  be 
found  in  [3]. 

5  Properties  of  the  conductivity-to-data  map  and 
preconditioning 

The  forward  map,  which  takes  conductivity  perturbations  67  to  voltage  data  w  through 
the  approximate  solution  of  (3)  will  be  denoted  by  E.  With  the  discretization  described 
previously,  this  map  is  an  m2  x  r»2  matrix,  operating  on  an  n2-vector  67  to  produce  an 
m2-vector  w.  To  construct  this  matrix,  we  use  the  Green’s  function  discussed  in  Section 
3  and  numerical  quadrature. 

With  this  notation,  the  inverse  problem  we  wish  to  solve  can  be  stated  as  a  system 
of  linear  equations  in  67 

E67  =  w.  (27) 

Here,  the  vector  w  is  the  measured  (or  synthetic)  data. 

Numerical  results  reported  in  [3],  and  our  own  experimentations  with  the  filtered 
Barber-Brown  backprojection  on  synthetic  data,  lead  us  to  believe  that  B  is  a  crude 
approximate  inverse  to  E.  Thus  we  are  led  to  consider  an  alternate  problem 

BE  67  =  Bw.  (28) 

How  well  one  can  solve  (27)  and  (28)  depends  on  the  properties  of  E  and  BE ,  and 
the  method  employed  for  the  solution.  We  are  ultimately  interested  in  solving  the  full 
nonlinear  problem  described  in  Section  2.  A  reasonable  approach  to  such  a  problem  is  to 
use  a  Gauss-Newton  method,  where  at  each  step,  we  need  to  solve  a  linear  system  very 
much  like  (27).  However,  this  step  is  necessary  only  to  obtain  an  update  towards  the  final 
solution.  Thus  in  the  early  stages  of  the  Gauss-Newton  method  it  is  often  sufficient  to 
solve  the  resulting  linear  system  only  approximately  (cf.  [7]).  With  this  in  mind,  we  rule 
out  direct  methods  for  the  solution  of  (27)  and  (28),  and  consider  only  iterative  methods. 
Of  particular  interest  are  conjugate  direction  algorithms  [12],  which  we  expect  will  yield 
a  good  approximate  solutions  to  (27)  and  (28)  in  a  small  number  of  iterations. 

Since  E  will  not  in  general  be  square,  it  is  natural  to  consider  the  normal  equation 

EtE  67  =  Etw  (29) 
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in  order  to  use  a  conjugate  direction  method.  We  could  in  principle  use  an  algorithm  of 
Hestenes  [13]  which  avoids  forming  the  normal  matrix  explicitly.  However,  regardless  of 
the  choice  of  the  algorithm  the  problem  solved  is  that  of  equation  (29).  The  conditioning 
of  this  solution  procedure  (convergence  rate  of  the  iterations  and  sensitivity  of  the  solution 
to  noise  in  the  data)  depends  on  the  eigenvalues  of  ET E. 


Consider  now  the  problem  (28).  BE  is  a  square  matrix  of  order  n2,  but  it  is  of  course 
not  symmetric,  which  leads  to  difficulties  with  many  conjugate  direction  methods.  The 
conjugate  residual  method  (CR)  is  a  special  example  of  a  conjugate  direction  method 
which  is  guaranteed  to  converge  for  any  square  system,  provided  the  matrix  has  definite 
symmetric  part.  As  will  be  seen  later,  the  symmetric  part  of  BE  is  not  exactly  definite, 
however  nearly  all  of  its  large  eigenvalues  have  same  sign  (positive).  In  practice  this 
should  insure  that  conjugate  residual  applied  to  BE  will  converge  for  a  wide  range  of 
initial  guesses.  From  the  point  of  view  of  computation,  this  method  is  advantageous 
because  it  works  on  the  matrix  BE  and  not  the  matrix  (BE)T BE.  The  eigenvalues  of 
BE  may  be  viewed  as  determining  the  conditioning  of  this  solution  procedure  [9], 


We  are  only  interested  in  the  action  of  E  and  BE  on  vectors  corresponding  to  images 
contained  in  the  unit  circle.  Therefore  we  remove  columns  of  E,  and  columns  and  rows 
of  BE  corresponding  to  pixels  lying  outside  the  unit  circle. 


For  m  =  n  =  16  we  have  256  pixels  of  which  only  20S  corresponds  to  locations  inside 
the  circle.  We  computed  the  eigenvalues  of  the  symmetric  part  of  the  matrix  BE  and 
found  that  of  the  the  moderate  to  large  size  eigenvalues,  all  except  for  one  are  positive 
(see  Figure  3a).  The  eigenvector  corresponding  to  the  negative  eigenvalue  is  displayed 
in  Figure  3b  and  seems  to  represent  a  constant  background  with  perturbations  near  the 
boundary.  The  difficulties  posed  to  the  conjugate  residual  method  by  the  presence  of  this 
"bad”  direction  can  be  avoided  if  one  restricts  6y  to  be  supported  sufficiently  far  away 
from  the  boundary. 


In  contrast,  the  symmetric  part  of  E  has  about  equally  many  large  positive  and 
large  negative  eigenvalues.  This  explains  why  even  in  the  square  case  we  cannot  apply 
conjugate  residual  successfully  to  E  alone. 


As  pointed  out  earlier,  we  believe  that  the  eigenvalues  of  ETE  and  BE  provide  in¬ 
formation  about  the  conditioning  of  our  solution  procedures  for  problems  (27)  and  (28) 
respectively.  The  eigenvalues  of  ETE  are  displayed  in  Figure  4.  The  fact  that  the  eigen¬ 
values  beyond  the  62th  are  all  less  than  10-3  reflects  the  illposedness  of  the  problem. 
The  largest  eigenvalue  is  1.958.  The  moduli  of  the  eigenvalues  of  BE  are  shown  in  Fig¬ 
ure  5  (BE,  being  nonsymmetric,  has  some  small  complex  eigenvalues  starting  with  the 
73rd).  We  find  that  the  conditioning  is  somewhat  better  here;  the  largest  modulus  of  any 
eigenvalue  is  1.080  and  eigenvalues  whose  moduli  are  less  than  10~3  begin  with  the  87th. 
The  tapering  off  of  the  eigenvalues  in  both  cases  indicates  that  some  resolution  loss  in 
reconstructing  6y  is  inevitable. 


Let  gj  denote  the  eigenvectors  of  ETE,  with  the  corresponding  eigenvalues  in  de¬ 
scending  order.  The  iterates  of  conjugate  direction  algorithms  in  practice  often  appear  to 
be  related  to  projections  of  the  true  solution  onto  the  span  of  ,j  for  increasing 

J .  Therefore  it  is  instructive  to  consider  the  projections  of  a  fixed  vector  onto  the  span 
of  }y— i, Let  hj  denote  the  eigenvectors  of  BE.  We  also  consider  the  projections 
of  the  same  fixed  vector  onto  the  span  of  { }y=i ,  ,j- 


We  take  a  profile  whose  entries  are  zero  except  for  one  corresponding  to  the  pixel 
whose  midpoint  is  (0.0625,  0.0625).  We  project  this  profile  onto  the  space  spanned  by 
the  {9j}j= i,  ,j  for  J  =  30  and  J  =  100.  The  results  are  shown  in  Figures  6a  and  6b. 


It  is  not  clear  if  any  iterative  method  will  get  as  far  as  J  =  100  because  the  100th 


S IBHHiK: 


eigenvalue  of  ET E  is  of  the  order  10"7.  In  some  sense  Figure  6b  gives  an  indication  of  the 
the  minimal  resolution  loss  to  be  expected.  Notice  that  with  J  =  30,  the  projected  image 
is  very  far  from  the  point  profile.  Indeed  the  Barber-Brown  method  produces  a  superior 
image  using  noiseless  synthetic  data  generated  from  the  point  profile  (compare  Figure 
6a  and  Figure  7).  There  are  two  eigenvalues  of  BE  with  nonzero  imaginary  part  among 
the  first  100  (the  conjugate  pair  corresponds  to  number  73  and  74).  The  projection  of 
any  real  vector  onto  <  100,  must  therefore  be  real  except  for  J  =  73. 

We  project  the  same  point  profile  as  before  onto  the  space  spanned  by  for 

J  =  30, 100.  In  both  cases  the  projected  image  is  superior  to  the  projected  image  using 
the  eigenvectors  of  ETE.  (compare  Figure  6  with  Figure  8). 

We  believe  that  this  numerical  study  indicates  that  the  iterates  of  CR  for  equation 
(28)  will  converge  faster  and  ultimately  get  closer  to  the  a  consistent  profile  than  the 
iterates  of  CR  when  applied  to  the  normal  equation  (29). 

To  separate  the  effects  of  the  filter  from  the  backprojection,  we  also  computed  the 
eigenvalues  of  the  symmetric  part  of  BE.  We  found  that  the  eigenvalues  are  similar  in 
structure  to  those  of  the  symmetric  part  of  BE  (only  one  large  negative  eigenvalue).  The 
presence  of  the  filter  does  however  increase  the  size  of  the  eigenvalues. 


0  20  40  60  80  100 

index 

Figure  4  Plot  of  the  distribution  of  the  eigenvalues  of  ET E.  The  largest  eigenvalue  is  1.958, 
and  the  eigenvalues  beyond  the  62nd  are  less  than  10-3. 


Figure  5  Plot  of  the  distribution  of  the  moduli  of  the  eigenvalues  of  BE.  The  largest  modulus 
is  1.080,  and  the  eigenvalues  beyond  the  87th  have  moduli  less  than  10-3. 
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Figure  6a  A  point  image  is  selected,  with  all  pixels  zero  except  at  location  (0.0625,  0.0625) 
which  has  value  1.  This  image  is  projected  onto  the  span  {g,,j  =  1, . . . ,  30}.  The  eigen¬ 
values  have  been  arranged  in  descending  order,  g}  is  the  >th  eigenvectors  of  ET  E.  The 
maximal  value  of  the  projection  is  0.0213,  the  minimal  value  is  -0.0170. 


Figure  6b  Projection  of  the  same  point  image  as  used  in  Figure  6a  onto  the  span  {g},j  = 
1, . . . ,  100}.  Maximal  value  =  0.2970,  minimal  value  =  -0.0522. 
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Figure  8a  Projection  of  the  same  point  image  as  used  in  Figure  6  onto  the  span  {hJtj  = 
1, . . . ,  30}.  Compare  to  Figure  6a.  Maximal  value  =  0.1033,  minimal  value  =  -0.0254. 


Figure  8b  Projection  of  the  same  point  image  as  used  in  Figure  6  onto  the  span  {hj,  j  = 
1, . . . ,  100}.  Compare  to  Figure  6b.  Maximal  value  =  0.5715,  minimal  value  =  -0.2105. 


6  An  iterative  refinement  for  the  reconstruction 


The  conjugate  residual  method  [8,  9]  solves  the  following  linear  equation: 

Ax  =  /, 

where  A  is  a  square,  but  possibly  nonsymmetric,  matrix  of  order  N  with  definite  sym¬ 
metric  part.  Suppose  the  initial  guess  is  xo  and  the  initial  residual  is  ro  =  f  —  Ax0.  The 
iterates  X{  are  ideally  obtained  by  minimizing 

\\A*i  ~  /II 

over  the  translated  Krylov  space 

*o  +  {ro,  Ar0,  •  ■  • ,  A,-1r0}. 

The  minimization  is  performed  by  consecutive  line  searches  in  directions  that  satisfy  the 
conjugacy  conditions 

(Apj,Api)  =  0,  for  j±i. 

In  practice,  for  a  nonsymmetric  matrix  A,  one  does  rarely  satisfy  all  the  conjugacy 
conditions  above.  The  version  of  the  algorithm  we  have  implemented  in  general  only 
guarantees  that  (Apj_i,  Api)  =  0.  This  version  is: 

(i)  choose  an  initial  guess  xq. 

(ii)  compute  the  initial  residual  r0  =  f  ~  A  x o- 

(iii)  set  search  direction  p0  =  ro- 

(iv)  for  i  =  0  step  1  until  convergence  do 

(a)  a,  =  (r,-,  Ap.j/fAp,,  Api) 

(b)  xi+l  =  n  +  a,p, 

(c)  ri+l  =  r,  -  aiApi 

(d)  bt  =  — (Arj+1,  Api)/{Api,  Api) 

(e)  P,+i  =  r,+1  +  biPi 

Although  this  version  of  the  conjugate  residual  method  does  not  guarantee  that  all  the 
conjugacy  conditions  (Apj,Api)  =  0,  j  ^  i,  are  satisfied,  and  therefore  in  general  does 
not  give  iterates  which  minimize  the  residual  over  the  relevant  translated  Krylov  spaces, 
this  version  is  known  to  converge  when  applied  to  matrices  with  definite  symmetric  part. 
We  refer  the  reader  to  [8,  9]  for  a  detailed  analysis. 

For  the  solution  of  (29),  we  simply  set  A  -  ET E,  and  as  the  intial  guess,  we  take  the 
filtered  Barber-Brown  backprojected  image. 

For  images  supported  near  the  boundary,  one  encounters  instabilities  by  a  direct 
application  of  CR  to  (28)  with  A  =  BE.  We  see  this  as  a  manisfestation  of  the  presence 
of  the  “bad”  search  direction  corresponding  to  the  large  negative  eigenvalue  (cf.  Figure 
3).  To  eliminate  this  difficulty,  we  set  to  zero  all  rows  and  columns  of  BE  corresponding 
to  pixels  outside  a  circle  with  a  priori  prescribed  radius  r  <  1.  In  matrix  notation,  this 
corresponds  to  post-  and  premultiplication  with  a  matrix  nr  obtained  from  the  identity 
matrix  by  setting  to  zero  the  columns  corresponding  to  pixels  outside  r  <  1.  In  our 
implementation  we  apply  CR  with  A  =  TlrBEUr,  for  an  appropriate  choice  of  r,  and 
we  take  the  filtered  Barber-Brown  backprojected  image  as  initial  guess.  If  the  image  6y 
is  not  supported  near  the  boundary,  one  may  use  r  =  1  without  any  difficulty.  If  r  is 
accidentally  chosen  too  small,  it  is  quite  easily  recognized  by  the  failure  of  the  residuals 
to  become  small. 
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7  Numerical  experiments 


We  conclude  this  paper  with  some  numerical  results.  In  Section  5,  we  made  a  prediction 
that  of  the  two  problems, 

EtE6j=Etw  (PI) 

BE  67  =  Bw,  ( /32) 

the  latter  is  more  well-behaved.  By  this  we  meant  that  the  iterates  of  the  CR  algorithm 
applied  to  equation  (P2)  should  converge  faster  and  ultimately  get  closer  to  a  consistent 
profile  than  iterates  of  the  same  CR  algorithm  applied  to  (PI). 

To  verify  our  prediction  we  show  the  results  of  computations  with  two  representative 
test  profiles.  In  both  cases  the  data  is  generated  by  a  multiplication  of  the  test  profile 
by  the  matrix  E,  i.e.,  by  modeling  perfect  dipoles  at  the  numbered  locations  in  Figure  2, 
and  solving  the  perturbational  equations  (3)  through  the  use  of  the  Green’s  function  and 
numerical  quadrature.  The  data  is  in  the  range  of  E  and  noiseless  (to  roundoff  errors). 

In  our  first  test  the  profile  used  is  shown  in  Figure  9a.  The  profile  represents  a  ring- 
shaped  high  conductivity  perturbation  with  a  ridge  across  (the  ridge  has  only  half  the 
strength  of  the  outer  ring).  The  ring  is  not  circular,  its  thickness  varies  and  it  is  off  center. 
In  Figure  10a  we  display  the  relative  problem  residual 

||E  tn  -  HI/IMI 

versus  number  of  iterations  of  the  CR  algorithm  for  both  equation  (PI)  and  (P2).  ||  J| 
denotes  the  Euclidian  norm.  It  is  clear  that  the  residuals  in  the  case  of  (P2)  are  smaller 
than  in  the  case  of  (PI)  (by  about  a  factor  of  1/2  at  the  30th  iterate).  This  is  also 
reflected  in  how  well  the  iterates  match  the  “correct”  profile  67.  In  Figure  10b  we  display 
the  relative  error 

ll*K-«7||/||*7ll, 

versus  number  of  iterations  of  the  CR  algorithm  for  both  equations  (PI)  and  (P2).  At 
the  30th  iterate  the  error  in  the  case  of  equation  (PI)  is  about  41%  where  as  in  the  case 
of  equation  (F2)  this  has  been  reduced  to  around  23%  (a  slightly  smaller  reduction  than 
for  the  residuals).  Most  impressive  to  observe  is  how  much  faster  the  CR  iterates  for  (P2) 
converge  during  the  first  6  steps  when  compared  to  those  for  (PI). 

For  further  comparison  we  examine  the  reconstructed  profiles  at  the  30th  iteration. 
For  reference  we  show  the  filtered  Barber-Brown  backprojection  (initial  guess  for  the 
iterative  schemes)  in  Figure  9b.  The  30th  iterates  for  (PI)  and  (P2)  are  show  in  Figures 
9c  and  9e  respectively.  Notice  that  the  ridge  is  recovered  in  the  case  of  (P2)  while  it  is 
still  not  visible  in  the  case  of  (PI).  In  fact  the  ridge  begins  to  be  visible  in  the  5th  CR 
iterate  of  (P2).  Figures  9d  and  9f  are  greylevel  plots  of  the  same  profiles  shown  in  Figures 
9c  and  9e  respectively.  The  ridge  is  clearly  visible  in  9f  whereas  it  is  not  seen  at  all  in 
9e.  In  the  computations  of  the  CR  iterates  for  (P2)  we  used  a  cutoff  radius  of  r  =  0.85 
to  avoid  instabilities,  as  described  in  the  last  part  of  Section  6. 

Our  second  test  involves  a  profile  in  the  form  of  two  rings,  one  twice  as  high  as  the 
other,  as  shown  in  Figure  11a.  We  display  the  relative  problem  residual  (Figure  12a)  and 
the  relative  error  (Figure  12b)  versus  the  iteration  number.  Again  the  CR  iterates  for 
(P2)  perform  significantly  better  than  those  for  (PI),  especially  in  the  first  few  steps. 
Figures  11b,  11c  and  lie  show  the  filtered  Barber-Brown  backprojection,  the  30th  iterate 
of  CR  applied  to  (PI)  and  the  30th  iterate  of  CR  applied  to  (P2).  Figures  lid  and  Ilf 
are  the  greylevel  plots  of  the  profiles  shown  in  Figures  11c  and  lie.  Notice  that  the  holes 
(the  areas  of  low  conductivity  inside  the  rings)  are  recovered  in  the  case  of  (P2)  (in  fact 
they  are  already  visible  in  the  7th  iterate),  while  they  are  still  invisible  in  the  case  of 
(PI).  Since  this  profile  is  supported  further  away  from  the  boundary  than  the  profile 


used  in  the  first  test  no  cutoff  is  necessary  in  the  computations  related  to  (P2)  (cf.  the 
end  of  Section  6) . 

These  two  tests,  and  several  others  we  have  performed,  indicate  that  B  acts  as  a  rea¬ 
sonably  good  preconditioner  for  the  original  problem.  It  is  fortunate  that  BE  is  positive 
definite  except  for  a  single  (controllable)  direction.  This  allows  us  to  use  CR  on  (P2). 
It  remains  to  be  seen  whether  BE  retains  the  same  property  when  E  corresponds  to  a 
background  conductivity  which  is  not  constant.  If  this  were  the  case,  then  we  can  con¬ 
struct  a  relatively  efficient  and  accurate  scheme  for  the  full  nonlinear  problem  (compared 
to  output  leastsquares). 

In  summary  we  conclude  that  the  CR  algorithm  applied  to  the  equation  BE6y  =  Bxv 
gives  an  iterative  method  which 

1)  allows  us  to  significantly  improve  upon  the  filtered  Barber-Brown  backprojection, 

2)  performs  better  than  the  conjugate  gradient  or  the  conjugate  residual  algorithm 
applied  to  the  normal  equation  ET  E6y  =  ETw 


(since  the  matrix  ET E  is  symmetric  and  positive  definite,  the  conjugate  gradient  and 
the  conjugate  residual  algorithms  will  behave  similarly).  We  are  hopeful  that  the  same 
result  may  be  obtained  for  much  more  general  domains. 


Figure  9a  The  test  profile 


Figure  9b  Reconstruction  using  the  filtered  Barber-Brown  backprojection. 
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Figure  9c  Reconstruction  based  on  (Pi)  at  the  30th  iterate  of  the  conjugate  residual  algorithm 
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Reconstruction  based  on  (P2)  at  the  30th  iterate  of  the  conjugate  residual  algo- 


Figure  Ilf  Gray  level  plot  of  Figure  lie. 
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Figure  12a  Reduction  in  the  relative  problem  residual  versus  number  of  iterations. 


Figure  12b  Relative  L 2  error  in  the  recovered  profile  versus  number  of  iterations. 
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